In many applications which are based on the SVD, arrival of new data requires SVD of the new matrix. Instead of computing from scratch, existing SVD can be updated.
The reader should be familiar with concepts of singular values and singular vectors, related perturbation theory, and algorithms.
The reader should be able to recognise applications where SVD updating can be sucessfully applied and apply it.
For more details see M. Gu and S. C. Eisenstat, A Stable and Fast Algorithm for Updating the Singular Value Decomposition and M. Brand, Fast low-rank modifications of the thin singular value decomposition and the references therein.
Let $A\in\mathbb{R}^{m\times n}$ with $m\geq n$ and $\mathop{\mathrm{rank}}(A)=n$, and let $A=U\Sigma V^T$ be its SVD. Let $a\in\mathbb{R}^{n}$ be a vector, and let $\tilde A=\begin{bmatrix} A \\ a^T\end{bmatrix}$. Then $$\begin{bmatrix} A \\ a^T\end{bmatrix} =\begin{bmatrix} U & \\ & 1 \end{bmatrix} \begin{bmatrix} \Sigma \\ a^TV \end{bmatrix} V^T. $$ Let $\begin{bmatrix} \Sigma \\ a^T V \end{bmatrix} = \bar U \bar \Sigma \bar V^T$ be the SVD of the half-arrowhead matrix. This SVD can be computed in $O(n^2)$ operations. Then $$\begin{bmatrix} A \\ a^T\end{bmatrix} = \begin{bmatrix} U & \\ & 1 \end{bmatrix} \bar U \bar\Sigma \bar V^T V^T \equiv \tilde U \bar \Sigma \tilde V^T $$ is the SVD of $\tilde A$.
Direct computation of $\tilde U$ and $\tilde V$ requires $O(mn^2)$ and $O(n^3)$ operations. However, these multiplications can be performed using Fast Multipole Method. This is not (yet) implemented in Julia and is "not for the timid" (quote by Steven G. Johnson).
If $m<n$ and $\mathop{\mathrm{rank}}(A)=n$, then $$ \begin{bmatrix} A \\ a^T\end{bmatrix} =\begin{bmatrix} U & \\ & 1 \end{bmatrix} \begin{bmatrix} \Sigma & 0 \\ a^T V & \beta\end{bmatrix} \begin{bmatrix} V^T \\ v^T \end{bmatrix}, $$ where $\beta=\sqrt{\|a\|_2^2-\|V^T a\|_2^2}$ and $v=(I-VV^T)a$. Notice that $V^Tv=0$ by construction. Let $\begin{bmatrix} \Sigma & 0 \\ a^T V & \beta\end{bmatrix} = \bar U \bar \Sigma \bar V^T$ be the SVD of the half-arrowhead matrix. Then $$\begin{bmatrix} A \\ a^T\end{bmatrix} = \begin{bmatrix} U & \\ & 1 \end{bmatrix} \bar U \bar\Sigma \bar V^T \begin{bmatrix} V^T \\ v^T \end{bmatrix} \equiv \tilde U \bar \Sigma \tilde V^T $$ is the SVD of $\tilde A$.
Adding a column $a$ to $A$ is equivalent to adding a row $a^T$ to $A^T$.
If $\mathop{\mathrm{rank}}(A)<\min\{m,n\}$ or if we are using SVD approximation of rank $r$, and if we want to keep the rank of the approximation (this is the common case in practice), then the formulas in Fact 1 hold approximately. More precisely, the updated rank $r$ approximation is not what we would get by computing the approximation of rank $r$ of the updated matrix, but is sufficient in many applications.
In [1]:
# pkg> add Arrowhead#master
using Arrowhead, LinearAlgebra
In [2]:
function mySVDaddrow(svdA::SVD,a::Vector)
# Create the transposed half-arrowhead
m,r,n=size(svdA.U,1),length(svdA.S),size(svdA.V,1)
T=typeof(a[1])
b=svdA.Vt*a
if m>=n || r<m
M=HalfArrow(svdA.S,b)
else
β=sqrt(norm(a)^2-norm(b)^2)
M=HalfArrow(svdA.S,[b;β])
end
# From Arrowhead package
svdM,info=svd(M)
# Return the updated SVD
if m>=n || r<m
return SVD([svdA.U zeros(T,m); zeros(T,1,r) one(T)]*svdM.V,
svdM.S, adjoint(svdA.V*svdM.U))
else
# Need one more row of svdA.V - v is an orthogonal projection
v=a-svdA.V*b
normalize!(v)
return SVD([svdA.U zeros(T,m); zeros(T,1,r) one(T)]*svdM.V,
svdM.S, adjoint([svdA.V v]*svdM.U))
end
end
Out[2]:
In [3]:
methods(SVD)
Out[3]:
In [4]:
import Random
Random.seed!(421)
A=rand(10,6)
a=rand(6)
Out[4]:
In [5]:
svdA=svd(A)
Out[5]:
In [6]:
norm(A*svdA.V-svdA.U*Diagonal(svdA.S))
Out[6]:
In [7]:
typeof(svdA)
Out[7]:
In [8]:
svdAa=mySVDaddrow(svdA,a)
Out[8]:
In [9]:
Aa=[A;transpose(a)]
println(size(Aa),size(svdAa.U),size(svdA.V))
[svdvals(Aa) svdAa.S]
Out[9]:
In [10]:
# Check the residual and orthogonality
norm(Aa*svdAa.V-svdAa.U*Diagonal(svdAa.S)),
norm(svdAa.U'*svdAa.U-I), norm(svdAa.Vt*svdAa.V-I)
Out[10]:
In [11]:
# Now flat matrix
Random.seed!(421)
A=rand(6,10)
a=rand(10)
svdA=svd(A);
In [12]:
A
Out[12]:
In [13]:
Aa
Out[13]:
In [14]:
svdAa=mySVDaddrow(svdA,a);
In [15]:
Aa=[A;transpose(a)]
println(size(Aa),size(svdAa.U),size(svdA.V))
[svdvals(Aa) svdAa.S]
Out[15]:
In [16]:
# Check the residual and orthogonality
norm(Aa*svdAa.V-svdAa.U*Diagonal(svdAa.S)),
norm(svdAa.U'*svdAa.U-I), norm(svdAa.Vt*svdAa.V-I)
Out[16]:
In [17]:
function mySVDaddcol(svdA::SVD,a::Vector)
X=mySVDaddrow(SVD(svdA.V,svdA.S,adjoint(svdA.U)),a)
SVD(X.V,X.S,adjoint(X.U))
end
Out[17]:
In [21]:
# Tall matrix
Random.seed!(897)
A=rand(10,6)
a=rand(10)
svdA=svd(A)
svdAa=mySVDaddcol(svdA,a);
In [22]:
# Check the residual and orthogonality
Aa=[A a]
norm(Aa*svdAa.V-svdAa.U*Diagonal(svdAa.S)),
norm(svdAa.U'*svdAa.U-I), norm(svdAa.Vt*svdAa.V-I)
Out[22]:
In [23]:
# Flat matrix
Random.seed!(332)
A=rand(6,10)
a=rand(6)
svdA=svd(A)
svdAa=mySVDaddcol(svdA,a);
In [24]:
# Check the residual and orthogonality
Aa=[A a]
norm(Aa*svdAa.V-svdAa.U*Diagonal(svdAa.S)),
norm(svdAa.U'*svdAa.U-I), norm(svdAa.Vt*svdAa.V-I)
Out[24]:
In [25]:
# Square matrix
A=rand(10,10)
a=rand(10)
svdA=svd(A);
In [26]:
svdAa=mySVDaddrow(svdA,a)
Aa=[A;transpose(a)]
norm(Aa*svdAa.V-svdAa.U*Diagonal(svdAa.S)),
norm(svdAa.U'*svdAa.U-I), norm(svdAa.Vt*svdAa.V-I)
Out[26]:
In [27]:
svdAa=mySVDaddcol(svdA,a)
Aa=[A a]
norm(Aa*svdAa.V-svdAa.U*Diagonal(svdAa.S)),
norm(svdAa.U'*svdAa.U-I), norm(svdAa.Vt*svdAa.V-I)
Out[27]:
In [28]:
# Adding row to a tall matrix
A=rand(10,6)
svdA=svd(A)
a=rand(6)
# Rank of the approximation
r=4
Out[28]:
In [29]:
svdAr=SVD(svdA.U[:,1:r], svdA.S[1:r],adjoint(svdA.V[:,1:r]));
In [30]:
# Eckart, Young, Mirsky
Ar=svdAr.U*Diagonal(svdAr.S)*svdAr.Vt
Δ=Ar-A
opnorm(Δ),svdvals(A)[5]
Out[30]:
In [31]:
svdAa=mySVDaddrow(svdAr,a);
In [32]:
Aa=[A; transpose(a)];
In [33]:
svdvals(Aa),svdvals([Ar;transpose(a)]),svdAa.S
Out[33]:
In [34]:
# Adding row to a flat matrix
A=rand(6,10)
svdA=svd(A)
a=rand(10)
# Rank of the approximation
r=4
Out[34]:
In [35]:
svdAr=SVD(svdA.U[:,1:r], svdA.S[1:r],adjoint(svdA.V[:,1:r]))
svdAa=mySVDaddrow(svdAr,a);
In [36]:
Ar=svdAr.U*Diagonal(svdAr.S)*svdAr.Vt
svdvals(Aa),svdvals([Ar;transpose(a)]),svdAa.S
Out[36]: